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ABSTRACT 

We show that in a hierarchical clustering model the low-order statistics of 
the density and the peculiar velocity fields can all be modelled semianalytically 
for a given cosmology and an initial density perturbation power spectrum 
P{k). We present such models for the two-point correlation function ^(r), the 
amplitude Q of the three-point correlation function, the mean pairwise peculiar 
velocity (f 12 the pairwise pecuhar velocity dispersion (fi2(^)), and the 
one-point peculiar velocity dispersion {vf). We test our models against results 
derived from N-body simulations. These models allow us to understand in 
detail how these statistics depend on P{k) and cosmological parameters. They 
can also help to interpret, and maybe correct for, sampling effects when these 
statistics are estimated from observations. The dependence of the small-scale 
pairwise peculiar velocity dispersion on rich clusters in the sample, for instance, 
can be studied quantitatively. There are also significant implications for the 
reconstruction of the cosmic density field from measurements in redshift space. 
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1. INTRODUCTION 

The large scale structure of the Universe is believed to have developed from small 
perturbations (usually assumed to be Gaussian) of the matter density field by gravitational 
instabilities. Under these assumptions the clustering pattern and velocity field observed 
today are determined by the initial conditions via the perturbation power spectrum [-P(fc)] 
and the cosmological parameters such as VLq, the cosmic density parameter. It is therefore 
possible to derive constraints on model parameters from the observed density and velocity 
distributions of galaxies. 
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There are two fundamental problems that must be addressed here: First, since it is 
unlikely that galaxies are completely Tinbiased testing particles of the matter density field, 
we need to understand any such bias in order to make meaningful comparisons between 
models and the observations of the galaxy distribution. Second, even if the observed galaxy 
distribution traces the matter distribution, we still need to understand how the observed 
distribution is related to the parameters that describe a cosmogonical model. The latter 
is by no means trivial, because the clustering pattern and velocity field observed today 
are nonlinear. N-body simulations are usually invoked to find a solution to this problem. 
However, as discussed comprehensively by Peebles (1980), a complete statistical description 
of the clustering process is also provided by the whole BBGKY hierarchy for the distribution 
function of mass particles. Indeed, adding some assumptions, one can solve the low order 
moments of the distribution function directly from the BBGKY equations. Although an 
incomplete description, such solutions are extremely useful for us to gain physical insights 
into the clustering process. These low order moments are also the most important ones, not 
only because they describe the most fundamental part of the clustering process but also 
because in practice they are the ones that can be measured from observations. 

In this paper, we show that semianalytical models can be constructed for all low order 
moments of the distribution function. These models allow such moments to be calculated 
directly from the initial density perturbation power spectrum for a given cosmological 
model, and thus provide us with a clear picture of how these moments are related to various 
model parameters. We present our models and test them against results from N-body 
simulations in Section 2. In Section 3 we demonstrate how our models can be used to help 
us to understand better some sampling effects when these moments are to be estimated 
from observations. Finally, in Section 4 we discuss the implications of our models to the 
reconstruction of real-space quantities from redshift distortions, and to the construction of 
some new statistics that can discriminate current models of structure formation. 



2. MODELS AND TEST BY N-BODY SIMULATIONS 

In this paper we will consider cosmogonies in which the Universe is dominated by cold 
dark matter (CDM). The cosmology is described by the cosmological matter density (Qq), 
the cosmological constant (Aq) and the Hubble constant (Hq = 100/ikms~^ Mpc~^). When 
Ao 7^ 0, we assume the universe to be flat so that r2o + Aq = 1. 

The initial power spectrum is [Bardeen et al. (1986)]: 



P{k) oc kT^{k), 



(la) 
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with 

T(k) = i< 

2Mq 



and 

Following Efstathiou et al. (1992), we have introduced a shape parameter, F = floh, for the 
power spectrum. (Note that the F defined here is not exactly the same as that of Efstathiou 
et al.). When F is treated as a free parameter, equation (1) can also be used to describe the 
power spectra in other structure formation models, such as the mixed dark matter (MDM) 
models (see e.g. Ma 1996). The rms mass fluctuation in top-hat windows with radius R, 
cr(i?), is deflned by 

a'{R)= r ^A\k)W\kR), (2a) 
Jo k 

where 

, 3(sina; — a; cosx) 
W{x) = ^ 3 (26) 

is the Fouriour transform of the top-hat window function, and 

A^k) = ^,k'P{k) . (3) 

For convenience, we will refer to A^(A;) as the power variance. We normalize P{k) by 
specifying erg = a{8 h~^Mpc). 



2.1. N-Body Simulations 

Before presenting our analytic approximations, we first give a brief summary of the 
N-body simulations to be used to test the models. 

We will use results derived from various P^M N-body simulations (see Jing et al. 1995 
for details on these simulations). Each simulation can be characterized by four model 
parameters (r^o; ^Oj T, cg), as discussed above, and three simulation parameters: the box 
size L (in h^^Mpc), the number of simulation particles Np and the effective force resolution 
e (in h^^Mpc). The model and simulation parameters for all simulations are listed in Table 
1. For each simulation we give it a name which is listed in the first column of the table and 
will be used throughout the paper to refer to the simulation. 
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2.2. The Two-Point Correlation Function 



The evolved two-point correlation function ^(r) is related to the evolved power variance 
A%(k) by 



Thus, in order to get ^(r) we need an expression for A^;(A;). Following the original 
argument of Hamilton et al. (1991), Peacock & Dodds (1994), Jain, Mo & White (1995), 
Padmanabhan et al. (1996) and Peacock & Dodds (1996, hereafter PD) have obtained 
fitting formulae which relate A^; to the initial density spectrum for a given cosmological 
model. The latest version of such a fitting formula is given in PD. PD have checked their 
fitting formula for Ae(A;) by a set of N-body simulations. Here wc use that formula and 
provide new tests for its predictions for ^(r) using results from (independent) N-body 
simulations. 

The evolved power variance Ag is assumed to be a function of its linear counterpart A: 

AiikE) ^ fmk,)], (5) 

where /c^ = [1 + A'^{kE)]~^^^kE- The functional form of / is approximated by (see PD for 
details) 



f{x) = X 



n{x) 



Q{x) 



where 



with 



Q{x) = A"'3[5(a)]3/^l^-^a;("-i/2)/3^ 

^ = 0.482 [l + n(A;L)/3]-°-9^^ 

5 = 0.226[l + n(A;L)/3]-^-^^^ 

y = 11.55[l + n(A;i)/3]-°-^'3 

a = 3.310[l + n(A;i)/3]-o-2^4 

/?- 0.862[l + n(A;L)/3]-°-'^^ 
, , , c? In P , , , , , 



(6) 

(7a) 
(76) 

(8a) 
(86) 
(8c) 

{8d) 
(8e) 

(8/) 
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The function g{a) in equation (7b) is the hnear growth factor at the cosmic time 
corresponding to the expansion factor a, and for given a-nd Aq, it is accurately described 
by (see Carroll, Press & Turner 1992): 

g{a) = [fi'^/^ - A + (1 + Q/2){1 + A/70)]"^ , (9a) 

where 

^ = ^(«) = FTTi ^° w . N ' (9^) 

^' a + Qo(l - a) + Ao(a3 - a) ^ ' 

A = A(a) = — — — ^ ^ , n r- (9c) 

^ ' a + Qo(l - a) + Ao(a3 - a) ^ ^ 

For later use we also define 

D{a) = ag{a). (9d) 

It is clear that A^; at any wavenumber is determined by A for a given cosmology, and 
we can just use equation (4) to obtain ^(r). In practice we start with the initial power 
variance A(/c) to obtain a table of [/cg, Ab(A;e)]. An interpolation is then used to carry out 
the integration in equation (4). 

The solid curves in Figure 1 show the model predictions for ^(r) in the SCDM models 
with (Ts = 1.24 and 0.62, and in the FLAT and OPEN models with (jg = 1. These 
predictions are compared to the results of N-body simulations shown by the symbols. For 
SCDM model with a% — 0.62, two different simulation box sizes are used to show the effects 
of both finite box size and finite simulation resolution on small scales. It is clear that the 
fitting formula works well for all cases. We have also obtained results for other models (i.e. 
FLAT models with Qq = 0.3 and 0.1, OPEN models with = 0.3 and 0.1, all having 
(78 = 1) and found a similar good agreement between model predictions and simulation 
results. 



2.3. The Mean Pairwise Peculiar Velocities 



From the pair conservation equation (Peebles 1980, §71), the ensemble (pair weighted) 
average of the pairwise pecuhar velocity (i'i2(r)) = ([v(x) — v(x + r)] • f) can be written as 



{vi2{r)) 



1 



diiy, a) 



Hr 3[1 + C(y,a)] 91na ' 

where r is the proper, and y the comoving, separation between the pairs; 



Ao(l - a"^) + Qo(a"^ - a"^) + a 



1/2 



(10a) 



(106) 
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is the value of Bubble's constant for an expansion factor a; 

a)^-J y'dy^y, a) = / -A|(fc, a)W{ky). (10c) 
y Jo Jo K 

Thus, in order to obtain (fi2(r)), we need to work out dAE{k,a)/da. Using the formulae 
presented in Subsection 2.2, it is straightforward to calculate the derivative dAE/da. In 
Appendix A, an explicit expression is given for this derivative. 

Figure 2 shows the comparison between the model prediction of {vi2{r)) and the 
simulation result. The agreement between the two is remarkably good for all cases except 
for B-SCDM0.62 where simulation result is significantly smaller than model prediction on 
r ^ h^^Mpc (the reason for this will be discussed below). This agreement shows again 
that the model of ^r) presented in §2.2 is also valid for describing the time evolution of 
i{r). Jain (1996) used the time evolution of the two-point correlation function derived 
directly from N-body simulations to solve for {vi2{r)) from the pair-conservation equation. 
The vahdity of our analytical model shows that {vi2{r)) can easily be calculated from the 
initial density spectrum via a single integration. The discrepancy between model prediction 
and simulation result for the case of B-SCDM0.62 is apparently due to the fact that this 
simulation is an intermediate output of SCDM1.24 and has evolved only by a factor of 4.5 
in the linear density growth factor since the initial simulation time. The numerical artifact 
of the initial density field generated by the Zel'dovich approximation should exist to some 
extent when the system is not sufficiently evolved. However, this artifact is expected to 
become smaller as the simulation evolves, as demonstrated by Baugh et al. (1995). It is 
therefore gratifying to see that the agreement between model prediction and simulation 
result is indeed better for SCDM1.24 than for B-SCDM0.62. 



2.4. Cosmic Energy Equation 



The (density weighted) mean square peculiar velocity of mass particles (vf) is related 
to the two-point correlation function by the cosmic energy equation: 



d 2/ 2\ A ^- ^^hia) 
—a^vf) = AnGpa^—^, 
da a m a 



where p is the mean density of the universe, and 



/2(a) = 



ydyiiy, a) 



dk A|(A;, a) 



Integrating equation (11) once, we have 

{vD = ln{a)H\a)a^l2{a) 



k 



0/2(0) Jo 



A;2 



l2{a)da 



(11) 



(12) 



(13) 
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In the linear case, hia) oc -D^(a) and 



1 - r D\a)da 

aD^{a) Jo ^ ' 



(14) 



We have found that equation (14) is a good approximation (to an error of < 10%) to 
equation (13) for all realistic cases. Thus, for a given cosmogonic model, we can easily 
obtain {v\). 

In Table 2 we compare the values of {v\) calculated from equation (14) to those 
obtained from N-body simulations. As one can see from equation (12), the contribution to 
{vD from the power at small wavcnumbcrs is not negligible, and the values of {vD derived 
from the simulations may be sensitive to the loss of power at k < 271 /L. To sec such an 
effect, we exclude from our analytic calculations all modes with k < 2tt/L. This is done 
by setting the lower hmit of the integration on the right hand side of equation (12) to be 
27r/L. The results are given in the left column under the column head '(v^)^/^ (model)'. 
For comparison, the values in brackets are obtained by assuming L — oo. As one can see, 
the effect of finite box size is indeed significant when the simulation box size is small. This 
effect is more important for the FLAT and OPEN models than for the SCDM model, 
because they have more power on large scales. Table 2 shows that the agreement between 
our model predictions and simulation results is reasonably good. 



2.5. Pairwise Peculiar Velocity Dispersion 

The relative velocity dispersion of particle pairs of separation r is defined as 
([v(x) — v(x + r)]^)^/^. In Fig. 3 we show (by symbols) the dispersion of the pairwise 
peculiar velocities projected along the separations of particle pairs ((t^i2(^))^'^^) in the 
N-body simulations. The main features of {vi2i''"))^^'^ are (a) monotonic rise at small r; (b) 
saturation at large r; (c) a maximum at medium r. In Section 2.5.1 we will discuss how 
these features can be explained by some physical arguments. Based on these arguments, 
we give (in Section 2.5.2) an empirical fitting formula for (f?2(^))^^^) so that it can be 
calculated easily for any given power spectrum. Readers who are mainly interested in using 
the model to calculate {vj2{r)y^'^ can skip Section 2.5.1 and go directly to Section 2.5.2. 
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2.5.1. Behaviour of the pairwise peculiar velocity dispersion 

At large separations where the correlated motion of particle pairs is negligible, the 
pairwise peculiar velocity dispersion is 

{vUr^oo))^l{vl). (15) 

Thus, the asymptotic value of {vl2{'>^)) is fixed at a constant for large pair separations. 

For very small separations, the main contribution to the pairwise velocity dispersion 
comes from particle pairs in virialized dark matter haloes (e.g. Marzke et al. 1994; Sheth 
1996; Shcth & Jain 1996). Assuming that dark haloes are spherically symmetric and that 
the velocities of particles in them are isotropic, we can write 

(^LW) = ^(^'(|x|) + ^^(|x + r|)), (16) 

where S{r) is the 3D velocity dispersion of particles at a radius r from a halo, and (■ ■ ■) 
denotes the pair- weighted average over all haloes. Thus, in order to obtain (fi2(^)) defined 
in equation (16), we need the density profile, as well as the mass function, of dark haloes. 

The results of N-body simulations (Hcrnquist 1990; Navarro, Frcnk & White 1996, 
hereafter NEW; Torman et al. 1996) suggest to write the density and velocity-dispersion 
profiles as 

n{r) = noF{r/rs), (17) 

S'{r) = S',{rs)G{r/r,), (18) 

where is a scale radius. For simplicity, wc assume that the velocity dispersion S{r) at a 
radius r is y^3/2 times the circular velocity at that radius, S^{r) = {3/2)GM{r)/r, where 
M(r) is the mass interior to r. Under such assumption G{x) is determined by F{x), 

G(x) = - r F{x)x^dx. (19) 
X Jo 

The total number of pairs of separation r in a halo can then be written as 

'P{r)dr = STT^lr^r^dr xldxiF{xi) J ^F{x2)d/3, (20) 

where 

xl^xi + {r/rsf - 2xi{r/rs)p. 
Similarly, the 5'^-weighted number of pairs in a halo is 

Vs{r)dr = 9>T:^nlrlSl{rs)r^dr xidxiF{xi) j ^ F{x2)[G{x^) + G{x2)]dl3. (21) 
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In carrying out the integrations in equations (20) and (21), we take F{r/rs) — when r is 
larger than the virial radius of the halo, Vy [which is defined as the radius of the mass shell 
that encloses the halo mass and settles at about half of its maximum expansion radius, see 
e.g. Lahav et al. (1991) for more detailed discussion]. 

We use the Press-Schechter formahsm (Press & Schechter 1974) to calculate the mass 
function of dark haloes. In this formalism, the comoving number density of dark haloes 
with mass in the range M ^ M -\- dM is 



/ 2\^/^ p Sc d\na(R) 



2(72 (i?) 



where M is the mass of the halo. M is related to the initial comoving radius R of the region 
from which the halo formed (measured in current units) by M = ^pR^, 5c is the threshold 
overdensity for collapsing. The predictions of equation (22) have been tested extensively by 
N-body simulations for various cosmogonies (e.g. Lacey & Coles 1994; Mo, Jing & White 
1996). 

With the above discussion, we can now write the pairwise velocity dispersion defined 
in equation (16) as: 

{Vi2[r)) - 3 ^oo N^M)V{r)dM ' ^ 

To complete the description, we need a model for Sc and the relation between the mass of 
a halo M and its virial radius r„. For the values of Sc in various cosmologies, we take the 
result summarized in Kochanek (1995; see also Bartclmann ct al. 1993). For the relation 
between r„ and M, we use the formulae given in Lahav et al. (1991). Thus, when Vg and 
the density profile F{r/rs) are given, we can easily obtain (fi2('")) defined in equation (23). 

In this paper we will use the model of NFW for the density profile of dark haloes: 

4^^ = ^[Ml + c) + (l+c)--ll-c3, (26) 

Pcrit -J 

where Pcrit is the current value of the critical density, m is the mass of a particle and c is 
the concentration factor, 

— ^200 /r,c\ 

c = , (26) 

with r2oo being the radius within which the average mass overdensity of the halo is 200. 
It should be pointed out that r2oo is in general different from the virial radius r^,, but for 
a given density profile the relation between and r2oo is fixed. The exact value of c for a 
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halo depends on its formation history and cosmology. The simulation result of NFW shows 
that c ~ 5-10 for typical haloes [i.e. those with Lagrangian radius it! ~ i?* where is 
defined by a[R^) = 1.69] in the standard CDM model. NFW also proposed a simple model 
for c based on halo formation time. The formation redshift Zf of a halo identified at redshift 
z < Zf with mass M is defined as the redshift by which half of its mass is in progenitors 
with mass exceeding f'M, where / < 1 is a constant. According to the formula given by 
Lacey & Cole (1993) based on the Press-Schechter formahsm, the halo formation time is 
then defined implicitly by 



erfc 



\, (27) 



where 5z is the threshold overdensity for collapsing at redshift z [e.g. = (1 + z)6c in an 
Einstein-de Sitter universe], a{M) is the rms of the linear power spectrum at redshift z — Q 
in top-hat windows enclosing mass M (see equation 2). NFW have used M200 (the mass 
enclosed in r2oo) instead of M in equation (27). We will follow their convention. NFW 
suggested that the characteristic overdensity of a halo identified at redshift z with mass M 
is related to its formation redshift Zf hy 

doiM, f, z) = Cimo [1 + Zf{M, f, z)]' , (28) 

where the normalization C{f) depends on /. We will take / = 0.01 as suggested by the 
N-body resuhs of NFW. In this case C{f) ^ 2 x lO''^ (NFW). Thus for a halo of given 
mass, one can obtain the concentration factor c from equations (25)- (28). In practice, we 
first solve Zf from equation (27) and insert the value of Zf into equation (28) to get 60, we 
then use this value of 6q in equation (25) to solve for c. As shown by NFW, the value of c 
does not depend sensitively on the choice of /, as long as / <C 0.1. It is necessary to point 
out that the arguments given in NFW are originally only for Einstein-de Sitter universe. 
Fortunately, these arguments are equally valid for low-Jl universes, as is shown recently by 
Navarro et al. (1996, in preparation). 

The dotted lines in Fig.3 show (t'i2(^")) given by equation (23) as a function of r. The 
pairwise velocity dispersion {vf2{r)) increases with r, because for large r the fraction of 
pairs from big haloes is larger. To examine the sensitivity of our results to the choice of 
the model for the concentration factor c, we also made calculations for the SCDM models 
by assuming c = 5 and 10. The results are shown in Figs. 3a and 3b as the dashed curves 
(how these curves are obtained will become clear in the following). A larger value of c gives 
higher amplitudes for (wi2(r)), because an increase of c enhances the velocity dispersion of 
a halo in the central region. 
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lt is clear that model (23) will break down for large r where the number of pairs within 
virialized haloes becomes smaller than other kinds of pairs (e.g. cross pairs between halo 
particles and field particles, and pairs among field particles). At some separation, the value 
of (^12) should decrease with r and approach the asymptotic value given by equation (15). 
Although the function {vl2{i^)) in the medium range of r is complicated, (because it depends 
not only on the amplitude but also on the shape of the three-point correlation functions, 
see below), it may still be possible to find some physically motivated models to understand 
the features observed in the simulations, given the asymptotic bahaviour of {vf2{''^)) on both 
small and large separations. 

We start with equation (72.1) in Peebles (1980): 

|(i + 0^" + ^i + 0^" + ^^(i + 0(^..^^./) 

^2^x" X" r^3^^^2Gpa [£xsC{l,2,3)^^0, (29) 

x"^ x"^ Jo J 

where v°' — {vu") is the a-component of the mean pairwise peculiar velocity (see equation 

10a), C(l,2,3) is the three point correlation function, m is the mass of a particle. The 

fourth term describes the mutual attraction of the particle pair and can be neglected when 

we consider dark matter particles. We write the velocity dispersion as 



6ap + [n - E] — -, (30) 



X 



where {vf) is the mean square peculiar velocity of particles (equation 11), 11 and E represent 
the effects of the correlated motions on the components of the dispersion parallel and 
perpendicular to the separation. Thus the (ID) velocity dispersion defined in the beginning 
of this subsection is {vl2{r)) = + If- Using the boundary condition that ^(r 00) = 
and n(r — * 00) = (see Peebles 1980, §72), we can obtain a formal solution for (t'i2(?^)) 
from equation (29): 



[1 + m Kir)) = M) + mHrfl - (Hrf 



d'^I dl d\nH dl 

+ 771 -TT. + 2- 



(?(lna)^ 9 In a 9 In a 5 In a 
dr Sn.H'^ f°° dr /■ 3 r ■ q 



+2/ -[l + e(r)](n-E) + ^— -Jd'q^ar,q,\r-q\), (31) 

where 



/ 

^0 



T^lMj^- (32) 

Note that r in equation (31) is the physical radius. Since models for ^, (vf) and have 
already been constructed, we need only to specify 11 — E and the integration of C in the last 
term of equation (31) to complete our model for (t'i2(r)). 
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For our demonstration, we will assume that the three point correlation function ( on 
small scales obeys the hierarchical form: 



(in, r2, rg) = Q [^(ri)e(r2) + ^(^2)^(^3) + ^(^3)^(^1)] , 



(33) 



where Q is a constant. It is clear from equation (31) that the pairwise peculiar velocity 
dispersion on small scales is dominated by the term of the three point correlation. In this 
case, we can write 



47r[l + e(r)] 



r J q-^ 



q|)+e(|r-q|)e(r)]. (34) 



Since ^ is known for a given model (see §2.2), the value of Q can then be determined by 
specifying {v\2ij)) at some small fiducial radius tq. On using our model for {v\2ij)) on small 
scales, the last term of equation (31) can therefore be fixed, as long as the hierarchical form 
(33) holds. On small scale, the hierarchical form is a reasonable approximation, however on 
large scales, Q is not a constant but depends on the size and shape of the triangle specified 
by (ri,r2,r3) (see e.g. Matsubara & Suto 1994; Jing & Borner 1996). Unfortunately, a 
theoretical model is not yet available to describe Q as a function of {ri,r2,r^) [|. In what 
follows, we assume equation (33) to hold on all scales. We will show that the error in 



'12 



V)) caused by such an assumption can be corrected by a simple model. 



To model 11 — S, we define a quantity which describes the anisotropy of the relative 
peculiar velocity dispersion: 

A[r) = 1 — / 2 / \\ — • (35a) 



It then follows that 



n 



{vUr)) 



A{r){vUr)). 



(356) 



On small scales where the velocity dispersion is isotropic, we have ^ = 0. On large scales 
where linear theory applies, we can write (see Peebles 1993, §21, note also the slight 
difference in the definitions of 11 and S): 



Mr) 



+ K2(r) 



(36a) 



r 3r3 3 ^ 



(366) 



^Note that the relevant quantity in our problem is the integration of ^ in the last term of equation (31), 
which is only a one-dimensional function of r and should therefore be much easier to model than itself. 
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where 



J3(r) = / r^^(r)dr = 




r"^ dk A'^{k) 
k [krY 



(sin kr — kr cos kr) ; 



JQ Jo 




kr\6 — {kr)'^ cos kr + 3 (kr)"^ — 2 sin kr 



K2{t) — I r^{r)dr — 




r-^ dk A'^{k) 
k {krY 



cos kr. 



In the above equations, A^(A;) is the dimensionless power spectrum in the Hnear regime. 
We will, however, insert the nonlinear power spectrum so that some nonlinear effects 
can be taken into account. Figure 5 compares A predicted by equations (35) and (36) 
(solid curves) to that derived from the simulations. The overall anisotropy in the velocity 
dispersion is quite small, and on separations r ^ 2 /i~^Mpc it is reasonably well described 
by our model. Since the (11 — E) term is expected to be important only on large scales, we 
will use equation (36) to estimate (11 — E) . 

The dot-dashed curves in Fig. 3 show the pairwise peculiar velocity dispersion 
(■^12 ('^))^''^ given by equation (31) with the assumptions discussed above. We have taken 
tq = 0.1 /i~^Mpc. As one can see from the figure, the result is not sensitive to the exact 
value of rg, as long as rg is small. Our model successfully gives a maximum in (i'i2(?^))^^^, as 
observed in the N-body simulations. The maximum occurs at a separation characteristic of 
the size of typical clusters, about 1-2 /i~^Mpc. This happens because most particles forming 
pairs of larger separations are in low-density environments where the velocity dispersion is 
low. The model overestimates the pairwise peculiar velocity dispersion on large scales when 
compared to the results derived from the N-body simulations (shown by the symbols). This 
overestimation comes mainly from our assumption that Q is a constant. Indeed, the value 
of Q may decrease substantially on large scales, as shown by N-body simulations (e.g. Jing 
& Borner 1996). The contribution to {v\2{t)) from the three-point correlation is therefore 
overestimated on large scales in our model. As discussed above, a reliable model for C, over 
different scales is still lacking at present time, and we have to invoke some simple arguments 
to correct the errors caused by the assumption of equation (33). 

As we have discussed above, for small separations r, most pairs are in virializcd haloes 
and equation (33) should be a reasonable description. The assumption of equation (33) fails 
for large r, because pairs outside dark haloes become important. As the simplest choice, the 
fraction of particles in dark haloes with radius r2oo > O-h''^ {<^h being a constant), fh{<^hi^), 
may be a good guess of the function that controls the change of (fi2('")) from its small 
separation value to its large separation value. According to the Press-Schechter formalism 
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(see equation [22]), fh{(y-hf) can be written in terms of the initial power spectrum as: 

5r 



fh{oihr) = erfc 



(37) 



where Rh is the hnear radius of a halo with r2oo = ct/i^- Thus when r is very small, all 
particles are in haloes so that fh = 1, while fh = when r oo. This behaviour of fh 
suggests the following empirical formula for 12 (?")): 

{vUr)) = [Maur)nvUr))s + {1 - [A(«.r)]^} x '^-{vj), (38) 

where we have changed our notation so that (i'i2(^))s denotes the {vf2{''^)) defined by 
equation (31). The values of ah and r are to be calibrated by results from N-body 
simulations. 

The thick sohd curve in Figure 3a is calculated according to (38) with the values 

= 0.2; r = 1. (39) 

It is clear that the model fits the simulation results for the SCDM model with as — 0.62 
reasonably well. The predictions of (vi2(?')) for other cosmogonical models are shown by 
the thick curves in the other panels of Fig. 3, using the same values of ah and r as given 
in equation (39). The figure shows that our model also works for these cosmogonies. For 
comparison, the two dashed curves in Figs. 3a and 3b show the same results obtained by 
assuming the concentration factor c = 5 (lower curve) and 10. 

The good agreement between the model predictions and the N-body results shows that 
the main features in (t'i2(r)) can indeed be explained by our simple physical considerations. 



2.5.2. Fitting formula for the pairwise peculiar velocity dispersion 

The model for {vl2{r)) given in Section 2.5.1, although having physical motivations, 
is not easy to implement. In this section we provide a much simpler fitting formula for 

Based on the N-body results and the discussion in Section 2.5.1 we make the following 
ansatz for {vl2{r)). 

{vUr)f' = ^'■'Hr,(l>{r/r,), (40a) 

where (f){x) is a universal function and Tc is a nonlinear scale. We choose Tc to be the virial 
radius of M* haloes, r*. For a given power spectrum, the hnear radius of M* haloes, Tq , is 
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given by a{rQ) — 1. The relation between r* and Tq in different cosmological models can be 
obtained by using the formulae given in Appendix B. We approximate the functional form 
of (p{x) by: 

. N ^ <Po.{l + Bx-^) + Ax-^ 
'^^"'^ 1 + \/[(7(a)]0-35[(](a)]0-2a:-("+^)' ^ ^ 

where </)oo = \J\{vlY/'^ / {HrSl^-^); A, B , V , a > {3 > and e > are constant. It follows 
that for r oo, {v'l2iT)) is forced to have the asymptotic value given by equation (15). 
For r — > 0, {v\2{^))^^'^ oc so that it increases with r as a power law. The dependence 
on the linear growth factor g{a) is included to take into account the fact that for a given 
power spectrum, haloes with the same linear radius may have different density profiles (or 
concentrations) in different cosmological models (see Section 2.5.1). The power of g{a) can 
be fixed by comparing the amplitudes of ('fi2('"))^^^ at small r for models (having the same 
initial power spectrum) with or without a cosmological constant. We found that a value 
of about 0.35 gives a reasonably good fit to our N-body results. The dependence on VL is 
included because for a given virial radius of a halo the ratio S/OP'^ (where S is the velocity 
dispersion of the halo) is larger in a lower-density universe. 

The sohd curves in Fig.4 show the predictions of our fitting formula with 

A = 58.67; S = -0.3770; y = 4.434; (40c) 

a = 2.25; (5 = 1.90; e - 0.15. (40ci) 

It is clear that the formula gives a reasonably good fit to all cosmogonical models examined 
in the present paper. At the moment, it is not clear whether the small discrepancy (typically 
^ 20 percent) between the model predictions and simulation results in some cases is due 
to the inaccuracy of our model, or due to the lower resolution of our simulations, or due to 
both of them. As one can see from Fig.4a, a lower resolution indeed reduces the pairwise 
peculiar velocity dispersion. Obviously, better simulations are needed to give a more 
accurate calibration of our formula. Given the uncertainties in our present simulations, we 
do not intend to obtain the best fit values for the parameters in equation (40b). 



2.6. The Hierarchical Amplitude of the Three-Point Correlation Function 

Given the model for 12 (?")), we can use equation (34) to calculate the hierarchical 
amplitude Q on small scales. In Figure 6 the solid curves show Q as a function of r 
predicted by our fitting formula (40). The amplitude Q decreases only slowly with r for 
r ~ 0.1/i~^Mpc, showing that equation (33) is vahd on small scales. The symbols in the 
figure show the values of Q estimated from N-body simulations by averaging over all triplets 
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with the smallest side equal to r and with the second smallest side less than 4r. Thus this 
quantity is not exactly the same as the one given by the model. However, since for small 
r the value of Q does not depend strongly on the shape of the triplets (see Jing 1996, in 
preparation), the comparison between the model and the simulation results may still be 
meaningful. The figure shows that our model prediction is generally in agreement with the 
simulation results. The values of Q at small separations (r ~ 1 /i"^Mpc) are typically 1-2, 
without a strong dependence on cosmogonies. These values of Q are only slightly larger 
than those obtained for galaxies (Q ~ 1 for optical galaxies, see e.g. Groth & Peebles 1977; 
Jing, Mo & Borner 1991). 

If the cosmic virial theorem is used under the assumptions that the two-point correlation 
function has a power-law form, ^{r) = A^r^'^, and that the three-point correlation function 
C has the hierarchical form given by equation (33), then the pairwise velocity dispersion on 
small scales can be written as: 

_ 2nGpQAor^-'rJi^) 
^"-^'^^-(7-l)(2-7)(4-7) ^ ^ 

(see Peebles 1980, Section 75). Here 

J(7)=/ dyy-\l + y''')M 
Jo 

with 

M^{1 + yf-^ [l - (2 - 7)y + y'] - |1 - y\^-^ [l + (2 - i)y + y\ . 

Thus, by using our model for {vf2{i^)) and by fitting ^{r) to a power law on small scales to 
get Ao, we can also get the value of Q from equation (41). The validity of equation (41) was 
checked by Suto (1993) using N-body simulations. There is, however, an uncertainty here. 
As one can see from Fig.l, ^(r) is not a good power-law even on small scales, the value of 
7 determined from ^ depends therefore on the range of r used in the fitting. Indeed, for 
r ~ 1 /i^^Mpc, the local value of 7 can sometimes be larger than 2 and equation (41) is 
ill-defined. Since it is unclear which range of r is most relevant for (t'i2(r)) on small scales, 
it does not seem to be a good idea to use 7 obtained from ^ in equation (41) to obtain Q. 
Fortunately, the value of 7 in equation (41) can also be determined from the dependence 
of {vf2{'>^)) on r. On small scales, {vf2{r)) increases with r approximately as a power law 
(Fig. 4). Thus the value of 7 given by fitting (fi2(?")) to a power law is always smaller than 
2, and equation (41) is well defined. The value of the correlation amplitude, Aq, can still 
be determined reasonably well by fitting ^(r) on small scales to a power law. As one can 
see from equation (41), the error in Q induced by Aq is only proportional to that in ^40. In 
Table 2, we show the values of Q predicted by equations (40) and (41) for ^(r) and {vf2{'r')) 
in the range r = 0.1 — 0.5 h~^Mpc, along with the values of 7 obtained by fitting (t'i2(r)) to 
a power law. The values of Q obtained in this way are similar to those shown in Fig. 6. 
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3. THE EFFECT OF REMOVING MASSIVE HALOES ON {vl^{r)) 

Davis & Peebles (1983) have calculated the pairwise velocity dispersion of galaxies for 
the CfAl redshift survey (Huchra et al. 1983). Their result, (t^^sl^))^^^ = 340 ±40kms-\ 
had been used for about 10 years to judge models of structure formation, and was a primary 
argument against the CDM model with f^o = 1 and with (Js ~ 1. Based on the same data 
set. Mo, Jing & Borner (1993, see also Zurek et al. 1994; Somerville, Davis & Primack 
1996) found that this statistic is sensitive to the presence (or absence) of galaxy clusters 
in a sample, and the constraint given by the observational results is uncertain. Recently, 
similar analyses have been performed on new, larger redshift surveys (Fisher et al. 1994; 
Guzzo et al. 1995; Marzke et al. 1995), but the problem remains: there is a large variation 
in (vi2('")) between different surveys, and even for the same survey analysed in different 
ways. From our model for {vl2{f)) (see Section 2.5.1), we see clearly why the observed small 
scale pairwise velocity dispersion is sensitive to the presence (or absence) of galaxy clusters. 
Indeed, from equations (20)- (23) we can write 

/-iV(M)MV3^M 
{Vu{r))oc ^^^^M^M^M ' ^^^^ 

where we have assumed that M oc and Sq oc M^/^. These assumptions are approximately 
correct for haloes with density profiles close to isothermal. Using equation (22) we see 
that while the number of pairs with small separations is dominated by dark haloes with 
M ~ (where is the mass at which the rms mass fluctuation o" = 1), the pairwise 
peculiar velocity dispersion on small scale can be significantly affected by massive ones 
with M > M*. The number density of such massive haloes is small and a large sample is 
therefore needed to have a fair samphng of their number density. Marzke et al. (1995) have 
discussed in considerable detail how such sampling can affect the statistics on (^12 ('"))• In 
this paper, we will not discuss this effect in detail. Instead, we show that from our model 
one can derive some new statistics which may be more robust against the sampling effect 
discussed above. 

Somerville, Primack & Nolthenius (1996) have suggested that studying (i'i2(r)) as 
clusters (dark haloes) with different internal velocity dispersions are removed leads to 
interesting information about the amount of power on cluster and subcluster scales. Since in 
practice the internal velocity dispersions of dark haloes are not easy to measure accurately 
for a large number of haloes, we suggest to study {vf2ii^)) as a function of the mean 
separation of the clusters that are removed from the analysis. 

From equation (22), we can write the total comoving number density of haloes with 
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mass exceeding M as 



N{> M) 



3 r"^ I 5c dhia 



exp 



51 1 dR 



(43) 



{2n)y^jR R^a{R)dlnR 



2(72 R ■ 



The mean separation of haloes with mass exceeding M is therefore d{M) = [A^(> M)]^^/^. 
As one can see from equation (43), for a given power spectrum, d determines R. Thus the 
effect of removing haloes with masses exceeding M in the model is to change the upper 
limits of the integrations in equation (23) from oo to M. This obviously has the effect of 
reducing (fi2(r)) on small scales. In Figure 7 we show {vl2{f))^^'^ at r = 0.1 h^^Mpc as a 
function of d for the SCDM models (F = 0.5) with erg = 0.62 and 1.24, and for three other 
cosmogonical models with F = 0.2 and as^o^ — 0.38. The use of the combination os^q^ 
is motivated by the fact that this quantity can be determined from observations based on 
large scale velocity fields (see e.g. Dekel 1994, Fisher et al. 1994) and on cluster abundance 
and clustering (e.g. White, Efstathiou & Frenk 1992; Mo, Jing & White 1996). The 
specific value 0.38 is chosen simply because the OPEN and FLAT models have this value. 
The figure shows that the value of {vl2{^))^^'^ does not reach its asymptotic value even for 
d ~ 50 /i~^Mpc (the typical value for the mean separation of rich clusters). Thus to have 
the small scale pairwise peculiar velocity dispersion fairly sampled, one needs a sample that 
contains many rich clusters. All galaxy redshift samples used so far to derive {vl2{'i^)) are 
not large enough to qualify for this purpose, and it is not surprising that the values of the 
small scale pairwise peculiar velocity dispersion of galaxies derived from available galaxy 
samples are uncertain. 

However, by inspecting Fig.7 we see that some useful statistics based on the pairwise 
velocity dispersion can still be derived from relatively small samples. For a given sample, 
we first identify from it clusters (groups) of different masses. We then remove successively 
the most massive clusters from the sample and analyze the small scale pairwise peculiar 
velocities as a function of the mean separation of clusters (d) that have been removed from 
the sample. Given a small sample, such a function can be determined reliably only for small 
values of d, because the mass function is well sampled only for relatively poor clusters. As 
one can see from Fig.7, the model prediction for this function is quite different for different 
cosmogonies. 

Another interesting point to note from Fig.7 is the difference in the results for the 
OPEN, FLAT and the (fio = l,r = 0.2, ag = 0.38) model. All these three models have the 
same shape of the power spectrum (F = 0.2) and the same value of ag^o^ (= 0.38). As we 
have discussed before, it is difficult to determine and erg separately from observations 
on large scale velocity fields or on cluster abundance. On the other hand, it is relatively 
easy to determine the combination crgJlo ^ from these observations. Fig.7 shows, however. 



-19- 



that one can hope to separate JIq and Ug by measuring the small scale peculiar velocity 
dispersion. The small scale pairwise peculiar velocity dispersion is lower in the model with 
f2o = 1 than in the OPEN and FLAT models with the same (TsJIq^, because clusters in this 
model formed late and thus are less concentrated. For the same reason, the small scale 
pairwise peculiar velocity dispersion is lower in the FLAT model than in the OPEN model. 
This line of argument is obviously worth further investigation. 

4. DISCUSSION 

The analytical models presented in Section 2 allow us to see clearly how the low order 
statistics of the density and the velocity distributions are determined by cosmological model 
and initial power spectrum. As a result we can use them to construct some statistical 
measures of the density and peculiar velocity fields to constrain models by observations. In 
this Section we will briefly discuss several possible applications of our results. Details of 
these applications will be described elsewhere. 

As a first application our results can be used to construct models for the reconstruction 
of cosmogonical parameters from measurements in rcdshift space. As discussed in Fisher 
et al. (1994, and references therein), in order to use the redeshift-space distortion on the 
two-point correlation function to derive cosmological parameters, it is necessary to model 
the distribution of the pairwise peculiar velocity. The distribution function is usually 
constructed from the low-order statistics of the pairwise peculiar velocity, which, by using 
our model, can all be calculated directly for any given cosmogony. Thus a liklihood function 
can be written for the pairwise peculiar velocities, which depends only on the cosmogonical 
parameters. Given a redshift sample of galaxies, one can then hope to constrain the 
cosmogonical parameters by maximizing the likelihood function. 

A second application of our models is to use them to examine what happens to some 
statistics, if the power spectrum is filtered in some manner. Such an investigation is 
important, because (1) real observations do not sample all parts of the spectrum equally 
well, and we need to understand how such sampling affects our statistics; and (2) we may 
sometimes want to filter the density field defiberately in order to get some insight into the 
underlying density field. In Section 3 we have already shown one such example related 
to the effect of removing massive clusters on the values of the pairwise peculiar velocity 
dispersions on small scales. Here we give another one. 

The mean square velocity (vf) given by the cosmic energy equation presented in 
Subsection 2.4 contains both large scale bulk motion and small scale random motion. It is 



sometimes desirable to separate these two kinds of motions by filtering the mass density 
spectrum. Prom equation (13), one can have a 'filtered version' of the cosmic energy 
equation: 



with W{k; x) being a filter with characteristic (comoving) radius of x. It is clear from 
equation (45) that if the high-frequency modes arc filtered out from the dimcnsionless power 
spectrum A%, then the velocity defined by equation (44) describes the motions induced by 
the low-frequency modes of the density perturbations. In linear case, this is just the mean 
square of the bulk motion of dark matter in the filter. On the other hand, if low frequency 
modes are filtered out, the velocity defined by equation (44) describes small scale relative 
motions. The statistics of the filtered velocity field can be derived either from velocity data 
(e.g. those based on Tully-Fisher relations) or from redshift-space distortions (see Miller, 
Davis & White 1996 for a discussion). Our model will then provide a clear picture how 
cosmogonical models are constrained by these measurements. 

It must be pointed out once again that in order to make any rigorous comparison 
between models and observations of galaxy distribution we need to understand how galaxy 
distribution is related to dark matter distribution. Much work needs to be done in this 
respect. However, for a given cosmogonical model, the density bias parameter, defined by 
the ratio between the two-point correlation function of galaxies and that of mass, and the 
velocity bias parameter, defined by e.g. the ratio between the pairwise peculiar velocity 
dispersion of galaxies and that of dark matter, are actually fixed. Our model can therefore 
be used to understand the effects of biased galaxy formation on the low-order moments of 
galaxy distribution with respect to those of dark matter. 

Finally, it is necessary to point out that the accuracy of our analytical models depends 
on the accuracy of the N-body simulations used to calibrate them. At present time, the 
discrepancy between our analytical formulae and simulation results is typically 20 percent. 
It is still unclear whether this small but significant discrepancy results from the inaccuracy 
of our analytical models, or from the limited resolution of N-body simulation, or from both 
of them. However, the method described in our paper is general. With improved N-body 
simulations, the fitting formulae for the evolved power spectrum and for the halo density 
profiles can all be improved. More accurate models for the low-order moments of dark 
matter distribution can then be easily obtained by using the procedure outlined in our 




(44) 



where 




(45) 



paper. 
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From equations (5)- (9) we see that A.%{k, a) depends on the expansion factor a through 
g{a) and x = A^(/ci), and we can write 



Princeton 



A. The Explicit Expression for dA\/da 




) 



k 




dg da dx \da ) 



{Al) 



This preprint was prepared with the AAS IM^jX macros v4.0. 
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The linear power spectrum evolves according to the linear theory so that x 
where Q{a) = [ag{a) /aog{ao)] and Xq = A'^^k^, oq) is the linear power spectrum 
expansion factor oq (e.g. at present time). Thus 



' dx^ 
da 



da 



dkr, \ da 



at a fiducial 



where /c^ = [1 + A|(A;, a)]~^/^/c. Inserting (A2) into (Al) we have: 



\ da , 



dldg^dldg^^- 
dg da dx da 



k 



dxo df 
din kL dx 



The derivatives on the right hand side of (A3) are: 



dg dg dD, dg dX 
da da dX da ' 

dfl _ _ fio(l - 1^0 - Ao + Sa^Ao) 
da [a + D.oi^ - a) + Xo{a^ - a)] 

dX _ a'^Xo[2a{l - Oq - Aq) + 31^o] 
da [a + Qq{1 - a) + Xo{a^ - a)] 

15r]^/Vl4 + 5/2-69A/28 



2 ' 



2 ' 



dg 
dn 

dg_ 
dX 



[1^4/7 _ A + (1 + n/2) (1 + A/70)] 

{bn/2) [!-(! + Vt/2) /70] 
[fi4/7 _ A + (1 + Vl/2) (1 + A/70)] 



2 ' 



2 ' 



d^_2^ 
da 



1 + 



rfln (yf' 



dl 

dg 



a \ cilna^ 
3/ Q{x) 



dx 



X 



1 + 



g 1 + Q(,t)' 
al3)Bx + an{x) (a - 1 /2) Q(a;) 



1 + 7^(a:) 1 + Q{x) 

The quantities in the above equations are defined in Subsection 2.2. 



(A3) 



{AA) 
{Ah) 

(AlO) 
{All) 
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B. The Relation between Linear and Virial Radii of Dark Haloes 

As shown by Lahav et al. (1991) the ratio between the virial radius and the 
turnaround radius of a dark halo can approximately be written as 

_ 1-77/2 . . 

where 

77 = 2AoQo'S-'(l + '^m)"', (B2) 

with Zm being the redshift at turnaround and S being the ratio between the halo and 
background densities at turnaround. The ratio between the linear radius (ro) and the virial 
radius (r^,) of a halo can then be written as 

The value of S depends on as well as on cosmological parameters. Viana & Liddle (1996) 
gave an accurate fitting formula for S as a function of fl = Q{zm)- 

where 

/(Q) = 0.76 - 0.25Q (open universe); (S5) 

/((]) = 0.73 - 0.23Q (flat universe). (S6) 

Thus, the relation between tq and is fixed for a given Zm- The value of Zm can be obtained 
by the assumption that the age of the universe at the turnaround of a halo {tm) is half the 
age of the universe at its collapse {tc): tm — tc/2. For a given redshift of halo collapse, 
Zci the time tc and the redshift z^ can be calculated from the redshift-time relation. This 
relation for an open universe is 

■ (1 + ^0^)^/^ Op , {i + nozY'' + {i-n,Yi' 

° {l-no){l+z) 2(l-fio)^/2 (l + fio^)l/2_(l_fi„)l/2' ^ ) 



and for a flat universe is 



2 , r[Ao(l + z)-3]V' + [Ao(l + ^)-^ + Jlo]'^' 
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Table 1. Simulation Parameters 



Model 


Realizations 


^^0 


Ao 


r 


(78 


L[/i-iMpc] 




e[/i-^Mpc] 


SCDM0.62 


3 


1.0 


0.0 


0.5 


0.62 


128 


100^ 


0.1 


B-SCDM0.62 


5 


1.0 


0.0 


0.5 


0.62 


300 


128^ 


0.2 


SCDM1.24 


5 


1.0 


0.0 


0.5 


1.24 


300 


128^ 


0.2 


FLAT 


5 


0.2 


0.8 


0.2 


1.00 


128 


643 


0.1 


OPEN 


5 


0.2 


0.0 


0.2 


1.00 


128 


643 


0.1 



Table 2. The Values of {vfy^^, Q and 7 



Model 




(model) 


{vfy^'^ (simulation) 


Q 


7 


SCDM0.62 


713 


[762] 


714 ± 22 


1.46 


1.74 


B-SCDM0.62 


742 


[762] 


694 ± 12 


1.46 


1.74 


SCDM1.24 


1464 


[1482] 


1380 ± 33 


1.55 


1.72 


FLAT 


521 


[595] 


524 ± 21 


1.78 


1.75 


OPEN 


571 


[642] 


556 ± 28 


1.42 


1.77 
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Fig. 1. — The two-point correlation functions of dark matter particles predicted by the 
analytical model (solid curves) compared to the results derived from N-body simulations 
(symbols). The dashed curves show the two-point correlation functions given by the linear 
power spectra. The error bars show the scatter among different realizations. The model 
parameters and the simulation box sizes are indicated in each panel. 
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Fig. 2. — The average pairwise peculiar velocity of dark matter particles predicted by the 
analytical model (solid curves) compared to the results derived from N-body simulations 
(symbols) . The error bars show the scatter among different reahzations 
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Fig. 3. — The pairwise peculiar velocity dispersion of dark matter particles predicted by the 
analytical model (curves) at the different stages of approximation (see Section 2.5.1). The 
dotted curves show the predictions of equation (23), while the dot-dashed and solid curves 
show the results given by equation (31) and equation (38), respectively. The two dashed 
curves (in a and b) show the predictions of equation (38) with the concentration factor c = 5 
(lower curve) and 10. The symbols show results derived from N-body simulations. Error 
bars are the scatter among different realizations 
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Fig. 4. — The pairwise peculiar velocity dispersion of dark matter particles given by the 
fitting formula (40) (curves), compared with the results derived from N-body simulations 
(symbols). Error bars are the scatter among different realizations 
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Fig. 5. — The velocity-dispersion anisotropy A of dark matter particles predicted by the 
analytical model (solid curves) compared to the results derived from N-body simulations. 
Each dotted curve shows the result of simulation in one particular realization. 
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Fig. 6. — The amplitude of the three-point correlation function of dark matter particles 
predicted by the analytical model (sohd curves) compared to the results derived from N- 
body simulations (symbols) . The error bars are the scatter among different realizations 
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Fig. 7. — The pairwise peculiar velocity dispersion of dark matter particles at separation 
r — 0.1 h~^Mpc as a function of d, the mean separation of haloes that are removed before the 
pairwise pecuhar velocity dispersion is calculated. Results are shown for the SCDM models 
(Qo = 1, r = 0.5) with (78 = 0.62 and 1.24, for FLAT and OPEN models (Qq = 0.2, F = 0.2) 
with CTg = 1 and for an additional model with = 1, F = 0.2 and erg = 0.38. The last three 
models have the same value of ctsOq^. 



